*          Program reg.test
*  To test the polynomial fitting subroutine 'reg'.
*
      character fnamein , fnameout , pname
      integer run
      real x(15),y(15),a(11)
      print,"Enter input data file name"
      read,fnamein
      print,"Output file name"
      read,fnameout
      open(unit=1,file=fnamein)

*     *** Read the run number and problem name

      read(1,1) run , pname
    1 format(i2,a16)

*     *** Read the number of data points and degree of
*     *** polynomial to be fitted.

      read(1,2) n , m
    2 format(2i3)

*     *** Read list of dependent and independent variables.

      read(1,3) (x(i),y(i),i=1,n)
    3 format(2f10.4)

      close(unit=1)
      call reg(pname,n,m,run,x,y,a,2,fnameout)
      end


      subroutine reg (pr , n , m , nplot , x , y , a , lfn , fname)
*
*  Purpose - To perform polynomial regression
*
*  Usage -
*            call reg (pr,pri,n,m,nplot,x,y,a,lfn,fname)
*
*  Arguments -
*
*      pr......Alphanumeric problem name.
*      n.......Number of observations.
*      m.......Highest degree polynomial specified.
*              If best fit polynomial is of degree
*              less than  m,  m  will be reset to
*              degree of best fit polynomial.
*      nplot...Option code for production of printer
*              plots of regression results.
*                  0  if plot not desired.
*                  1  if plot desired.
*      x.......vector of independent variables.
*      y.......vector of dependent variables.
*      a.......vector of regression coefficients.
*              There will be  m + 1  entries in  a .
*      lfn.....Logical unit number for output file.
*      fname...Output file name
*
*  Subroutines required - pplot , orthls , coeffs , fity
*
*  Remarks -
*      (1) The largest degree polynomial which may be fitted is 9.
*      (2) The number of data points must be greater than the
*          largest degree polynomial fitted.
*      (3) A maximum of 50 observations is assumed.
*
*--------------------------------------------------------------------
*
      character pr , fname
      real x(n) , y(n) , a(11)
      real w(50) , t1(50) , t2(50) , t3(50) , plt (150) , z(50)
      real c(10) , alpha(10) , beta(10) , b(10)
*
*--------------------------------------------------------------------
*
       open(unit=lfn,file=fname,access="sequential",recl=80)
*
     3 format("1Polynomial Regression.....",a16/)
     4 format(" Number of Observations",i6//)
     5 format(" Polynomial Regression of Degree",i3/1x,
&              "***********************************"/)
     6 format("   Coefficients"/3(1x,4e15.5/))
     8 format("0","      --- Analysis of Variance ---"/)
     9 format(" ",5x,"Source",8x,"DF",8x,"SS",11x,"MS",12x,"F"/)
    10 format(" ",2x,"Regression...",i6,3f13.3)
    11 format(3x,"Error........",i6,2f13.3)
    12 format(3x,"Total........",i6,f13.3/)
    13 format("   No Improvement")
    14 format(" "//27x,"Table of Residuals"//" Obs. No.",5x,
&             "X Value",6x,"Y Value",4x,"Y Estimate",4x,
&             "Residual"/)
    15 format(1x,i6,4f13.4)
    16 format(" ",11x,"Improvement in SS",f16.3///)
*
*--------------------------------------------------------------------
*
*     *** Print problem name and no. of obs., check order
*
      write(lfn,3) pr
      write(lfn,4) n
      if (m .gt. 9) m = 9
*
*     *** Find mean and total sum of squares for dependent variable
*
      ybar = 0.0
      ssat = 0.0
      do 100 i = 1 , n
         ybar = ybar + y(i)
         ssat = ssat + y(i) * y(i)
  100 continue
      fn = n
      ybar = ybar / fn
      ssat = ssat - fn * ybar * ybar
      nt = n - 1
*
*     *** Generate orthogonal polynomials and coefficients
*
      call orthls(x,y,w,n,0,0,c,alpha,beta,m,m+1,t1,t2,t3,nd1)
*
*     *** Check fit for successive degrees, and find sum of squares
*         attributable to error.  Compute -
*           ssar  - S of S attributable to regression
*           smar  - Mean square for regression
*           smae  - Mean square for error
*           fval  - F-Ratio
*           sumip - Improvement in SS
*           nr    - D of F for regression
*           ne    - D of F for error
*           nt    - D of F Total
*
      sum = 0.0
      do 200 k = 1 , m
         k1 = k
         call fity(x,n,0,c,alpha,beta,k,k+1,plt,t1,t2,iend3)
         ssae = 0.0
         do 120 i = 1 , n
            temp = y(i) - plt(i)
            ssae = ssae + temp * temp
  120    continue
         nr = k
         ne = nt - nr
         ssar = ssat - ssae
         smar = ssar / float(nr)
         smae = ssae / float(ne)
         fval = smar / smae
         sumip = ssar - sum
         sum = ssar
         write(lfn,5) k
         if (sumip) 140 , 140 , 150
  140       write(lfn,13)
            go to 210
  150    call coefs(0,c,alpha,beta,k,k+1,a,t1,t2,t3,ind2)
         k2 = k + 1
         write(lfn,6) (a(l),l=1,k2)
         write(lfn,8)
         write(lfn,9)
         write(lfn,10) nr , ssar , smar , fval
         write(lfn,11) ne , ssae , smae
         write(lfn,12) nt , ssat
         write(lfn,16) sumip
  200 continue
*
*     *** Plot or not
*
  210 if (nplot) 400 , 400 , 220
*
  220 call fity(x,n,0,c,alpha,beta,k2-1,k2,z,t1,t2,iend3)
      do 230 i = 1 , n
         plt(i) = x(i)
         i1 = n + i
         plt(i1) = y(i)
         i1 = i1 + n
         plt(i1) = z(i)
  230 continue
*
*     *** Print table of residuals
*
      write(lfn,3) pr
      k = k2 - 1
      write(lfn,5) k
      write(lfn,14)
      do 250 i = 1 , n
         resid = y(i) - z(i)
         write(lfn,15) i,x(i),y(i),z(i),resid
  250 continue
      call pplot(nplot,plt,n,3,3*n,n,1,68,lfn)
  400 m = k
      close(unit=lfn)
      return
      end


      subroutine orthls (x,y,w,n,l,j,c,alpha,beta,k,k1,t1,t2,t3,ind1)
*--------------------------------------------------------------------
*     This subroutine computes the coefficients of the polynomial
*     equation of degree k and the alpha and beta parameters.
*--------------------------------------------------------------------
      real x(n),y(n),w(n),c(k1),alpha(k1),beta(k1),t1(n),t2(n),t3(n)
*--------------------------------------------------------------------
*     Program initialization.
*--------------------------------------------------------------------
      kj1 = k - j + 1
      if (kj1 .le. 0) go to 16
      sum = 0.0
      if (l .eq. 1) go to 3
         do 2 i = 1 , n
            t3(i) = x(i)
            if (j .gt. 0) go to 1
               sum = sum + 1.0
               go to 2
    1       sum = sum + x(i)**(2*j)
    2    w(i) = 1.0
         go to 7
    3 do 6 i = 1 , n
         t3(i) = x(i)
         if (j .gt. 0) go to 4
            sum = sum + w(i)
            go to 5
    4    sum = sum + w(i) * x(i)**(2*j)
    5    x(i) = w(i) * x(i)
         y(i) = w(i) * y(i)
    6 continue
    7 b = 0.0
      r0 = sum
      do 9 i = 1 , n
         if (j .gt. 0) go to 8
            t2(i) = 1.0
            go to 9
    8    t2(i) = t3(i)**j
    9 t1(i) = 0.0
*--------------------------------------------------------------------
*     Begin computation.
*--------------------------------------------------------------------
      ii = 1
   10 s = 0.0
      do 11 i = 1 , n
   11 s = s + y(i)*t2(i)
*--------------------------------------------------------------------
*     Computation of a coefficient in the polynomial equation.
*--------------------------------------------------------------------
      c(ii) = s / r0
      if (ii .ge. kj1) go to 15
*--------------------------------------------------------------------
*     Computation of an alpha for the polynomial equation.
*--------------------------------------------------------------------
      sumxps = 0.0
      do 12 i = 1 , n
   12 sumxps = sumxps + x(i) * t2(i) * t2(i)
      alpha(ii) = sumxps / r0
*--------------------------------------------------------------------
*     Computation of a new polynomial.
*--------------------------------------------------------------------
      do 13 i = 1 , n
         temp  = t2(i)
         t2(i) = (t3(i) - alpha(ii)) * t2(i) - b * t1(i)
         t1(i) = temp
   13 continue
*--------------------------------------------------------------------
*     Computation of a beta for the polynomial.
*--------------------------------------------------------------------
      r = 0.0
      do 14 i = 1 , n
   14 r = r + w(i) * t2(i) * t2(i)
      beta(ii) = r / r0
      r0 = r
      b  = beta(ii)
      ii = ii + 1
      go to 10
*--------------------------------------------------------------------
*     Successful return.
*--------------------------------------------------------------------
   15 ind1 = 1
      return
*--------------------------------------------------------------------
*     Error return.  Set all c coefficients, alpha and beta to zero.
*--------------------------------------------------------------------
   16 do 17 ii = 1 , k
         c(ii) = 0.0
         alpha(ii) = 0.0
         beta(ii)  = 0.0
   17 continue
      c(k+1) = 0.0
      ind1 = -1
      return
      end


      subroutine coefs(j,c,alpha,beta,kc,k1,a,t1,t2,t3,ind2)
*-------------------------------------------------------------------
*     This subroutine computes the  a  coefficients for a
*     polynomial of degree  kc  where  kc  is less than or
*     equal fo  k .
*-------------------------------------------------------------------
      real c(k1),alpha(k1),beta(k1),a(k1),t1(k1),t2(k1),t3(k1)
*-------------------------------------------------------------------
*     Program initialization.
*-------------------------------------------------------------------
      kcj1 = kc - j + 1
      if (kcj1 .le. 0) go to 9
      b = 0.0
      do 1 nn = 1 , kcj1
         a(nn)  = c(nn)
         t1(nn) = 0.0
         t2(nn) = 0.0
         t3(nn) = 0.0
    1 continue
      if (kc .le. j) go to 5
      ii = 2
*-------------------------------------------------------------------
*     Begin computation.
*-------------------------------------------------------------------
    2 t2(ii) = 1.0
      do 3 nn = 2 , ii
         t3(nn)  = t2(nn-1) - t2(nn) * alpha(ii-1) - b * t1(nn)
         a(nn-1) = a(nn-1) + c(ii) * t3(nn)
    3 continue
      if (ii .ge. kcj1) go to 5
*-------------------------------------------------------------------
*     Reset the vectors for the next coefficient.
*-------------------------------------------------------------------
      do 4 nn = 1 , ii
         t1(nn) = t2(nn)
         t2(nn) = t3(nn)
    4 continue
      b  = beta(ii-1)
      ii = ii + 1
      go to 2
    5 if (j .le. 0) go to 8
*-------------------------------------------------------------------
*     Arrange coefficients properly if  j  is non zero.
*-------------------------------------------------------------------
      do 6 nn = 1 , kcj1
         n1 = kcj1 - nn + 1
         n2 = n1 + j
         a(n2) = a(n1)
    6 continue
      do 7 nn = 1 , j
    7 a(nn) = 0.0
*-------------------------------------------------------------------
*     Successful return
*-------------------------------------------------------------------
    8 ind2 = 2
      return
*-------------------------------------------------------------------
*     Error return, set all  a  coefficients to zero
*-------------------------------------------------------------------
    9 do 10 nn = 1 , kc
   10 a(nn) = 0.0
      a(kc+1) = 0.0
      ind2 = -2
      return
      end


      subroutine fity(xf,m,j,c,alpha,beta,kf,kf1,yf,t1,t2,ind3)
*---------------------------------------------------------------------
*     This subroutine computes the fitted values for a
*     given set of arguments with a polynomial of degree
*     kf where kf is less than or equal to k.
*---------------------------------------------------------------------
      real xf(m),c(kf1),alpha(kf1),beta(kf1),yf(m),t1(m),t2(m)
*---------------------------------------------------------------------
*     Program initialization.
*---------------------------------------------------------------------
      kfj1 = kf - j + 1
      if (kfj1 .le. 0) go to 7
      b = 0.0
      do 2 i = 1 , m
         yf(i) = 0.0
         if (j .gt. 0) go to 1
            t2(i) = 1.0
            go to 2
    1    t2(i) = xf(i)**j
    2 t1(i) = 0.0
*---------------------------------------------------------------------
*     Begin computation.
*---------------------------------------------------------------------
      ii = 1
    3 do 4 i = 1 , m
    4 yf(i) = yf(i) + c(ii)*t2(i)
      if (ii .ge. kfj1) go to 6
*---------------------------------------------------------------------
*     Computation of new polynomial.
*---------------------------------------------------------------------
      do 5 i = 1 , m
         temp = t2(i)
         t2(i) = (xf(i) - alpha(ii)) * t2(i) - b * t1(i)
         t1(i) = temp
    5 continue
      b  = beta(ii)
      ii = ii + 1
      go to 3
*---------------------------------------------------------------------
*     successful return
*---------------------------------------------------------------------
    6 ind3 = 3
      return
*---------------------------------------------------------------------
*     Error return. Set all fitted values equal to zero.
*---------------------------------------------------------------------
    7 do 8 i = 1 , m
    8 yf(i) = 0.0
      ind3 = -3
      return
      end


      subroutine pplot(n0,a,n,m,nm,nl,ns,nchar,lfn)
*---------------------------------------------------------------------
*       Purpose -
*            Plot several cross-variables versus a base
*            variable in strip chart form.  This is not
*            a tidy program and may have bugs in the y-axis
*            (looking sideways at the plot) labels.  Main use
*            was for quick looks at long strip chart like
*            arrays of data.  This program was evidently
*            taken from an old IBM SSP package.
*
*       Usage -
*            call pplot(n0,a,n,m,nm,nl,ns,nchar,lfn)
*
*       Arguments -
*            n0      - Chart number (3 digits maximum)
*            a       - Matrix of data to be plotted.  First
*                      column represents base variables and
*                      successive columns are the cross-variables
*                      (maximum is 9)
*            n       - Number of rows in matrix a
*            m       - Number of columns in matrix a
*                      (equal to the total number of variables)
*                      Maximum is 10
*            nm      - The total number of elements in a (usually n*m)
*            nl      - Number of lines in the plot.  If 0 is
*                      specified, 50 lines are used
*            ns      - Code for sorting the base variable data
*                      in ascending order
*                        0  sorting is not necessary (already in
*                           ascending order)
*                        1  sorting is necessary
*            nchar   - Number of print positions available
*            lfn     - Logical file number for printed output
*---------------------------------------------------------------------
      real      ypr(11) , a(nm)
      character blank   , ang(9) , out(121) , point

    1 format("1",20x,"CHART  ",i3,/)
    2 format(" ",f8.3,3x,111a1)
    3 format(" ")
    7 format(" ",11x,".",10(9x,a1))
    8 format(" ",5x,11f10.4)

* Initialize Parameters

      blank = " "
      point = "."
      do 5 i = 1 , 9
        ang(i) = cnvi2c(i)
    5 continue
      nchars = nchar/10
      nchars = 10*nchars
      chars  = nchars
      fy     = chars/10.0
      ny     = fy + 1.0
      nt     = nchars
      nll    = nl

* Perform sort or transposition as necessary

      if (ns) 18 , 18 , 10
   10 do 14 i = 1 , n
        do 13 j = i , n
          if ( a(i) - a(j)) 13 , 13 , 11
   11     l = i - n
          ll = j - n
          do 12 k = 1 , m
            l     = l + n
            ll    = ll + n
            f     = a(l)
            a(l)  = a(ll)
            a(ll) = f
   12     continue
   13   continue
   14 continue
      go to 18

* Test number of lines

   18 if (nll) 20 , 19 , 20
   19 nll = 50
   20 write(lfn,1) n0

* Find scales for base and cross variables

      xscal = (a(n)-a(1))/(float(nll-1))
      m1 = n + 1
      m2 = m*n
      ymin = a(m1)
      ymax = ymin
      do 40 j = m1 , m2
        if (a(j) - ymin) 26 , 28 , 28
   26     ymin = a(j)
   28   if (a(j) - ymax) 40 , 40 , 30
   30     ymax = a(j)
   40 continue
      yscal  = (ymax-ymin)/chars
      min    = ymin/(yscal*fy)
      ynomin = fy*yscal*float(min)
      if ( ymin .lt. ynomin )  ynomin = ynomin - fy*yscal
      ymin   = ynomin
      ymax   = ymin + chars*yscal

* Find base variable print position

      xb   = a(1)
      half = xscal/2.0
      l    = 1
      my   = m - 1
      do 80 i = 1 , nll
        f = i - 1
        xpr = xb + f*xscal
   48   if ( a(l) .ge. xpr-half ) go to 49
*         base variable too small
          l = l + 1
          go to 48
   49   if ( a(l) .gt. xpr+half ) go to 70

* Find cross-variable

        do 55 ix = 1 , nt
   55   out(ix) = blank
        do 60 j = 1 , my
          ll = l + j*n
          jp = ((a(ll)-ymin)/yscal) + 1.5
          out(jp) = ang(j)
   60   continue

* Print line and clear or skip
      l    = 1
        write(lfn,2) xpr , ( out(iz),iz=1,nt)
        l = l + 1
        go to 80
   70   write(lfn,3)
   80 continue

* Print cross-variable numbers

      write(lfn,7) (point,ip=1,ny)
      ypr(1) = ymin
      do 90 kn = 1 , ny-2
   90 ypr(kn+1) = ypr(kn) + yscal*chars/(ny-1)
      ypr(ny) = ymax
      write(lfn,8) (ypr(ip),ip=1,ny)

      return
      end
